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Asymptotically Optimal Algorithms for Pickup 
and Delivery Problems with Application to 
Large-Scale Transportation Systems 

Kyle Treleaven, Marco Pavone, Emilio Frazzoli 
Abstract 

Pickup and delivery problems (PDPs), in which objects or people have to be transported between specific locations, 
are among the most common combinatorial problems in real-world logistical operations. A widely-encountered type 
of PDP is the Stacker Crane Problem (SCP), where each commodity/customer is associated with a pickup location 
and a delivery location, and the objective is to find a minimum-length tour visiting all locations with the constraint 
that each pickup location and its associated deUvery location are visited in consecutive order. The SCP is NP-Hard 
and the best known approximation algorithm only provides a 9/5 approximation ratio. The objective of this paper 
is threefold. First, by embedding the problem within a stochastic framework, we present a novel algorithm for the 
SCP that: (i) is asymptotically optimal, i.e., it produces, almost surely, a solution approaching the optimal one as the 
number of pickup s/deliveries goes to infinity; and (ii) has computational complexity 0{n^^^), where n is the number 
of pickup/delivery pairs and £ is an arbitrarily small positive constant. Second, we asymptotically characterize the 
length of the optimal SCP tour. Finally, we study a dynamic version of the SCP, whereby pickup and delivery requests 
arrive according to a Poisson process, and which serves as a model for large-scale demand-responsive transport (DRT) 
systems. For such a dynamic counterpart of the SCP, we derive a necessary and sufficient condition for the existence 
of stable vehicle routing policies, which depends only on the workspace geometry, the stochastic distributions of 
pickup and delivery points, the arrival rate of requests, and the number of vehicles. Our results leverage a novel 
connection between the EucUdean Bipartite Matching Problem and the theory of random permutations, and, for the 
dynamic setting, exhibit novel features that are absent in traditional spatially-distributed queueing systems. 

I. Introduction 

Pickup and delivery problems (PDPs) constitute an important class of vehicle routing problems in which objects 
or people have to be transported between locations in a physical environment. These problems arise in many 
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contexts such as logistics, transportation systems, and robotics, among others. Broadly speaking, PDPs can be 
divided into three classes (T]: 1) Many-to-many PDPs, characterized by several origins and destinations for each 
commodity/customer; 2) one-to-many-to-one PDPs, where commodities are initially available at a depot and are 
destined to customers' sites, and commodities available at customers' sites are destined to the depot (this is the 
typical case for the collection of empty cans and bottles); and 3) one-to-one PDPs, where each commodity/customer 
has a given origin and a given destination. 

When one adds capacity constraints to transportation vehicles and disallows transshipments, the one-to-one PDP 
is commonly referred to as the Stacker Crane Problem (SCP). The SCP is a route optimization problem at the core 
of several transportation systems, including demand-responsive transport (DRT) systems, where users formulate 
requests for transportation from a pickup point to a delivery point Q, |[3|. Despite its importance, current algorithms 
for its solution are either of exponential complexity or come with quite poor guarantees on their performance; 
furthermore, most of the literature on the SCP does not consider the dynamic setting where pickups/deliveries are 
revealed sequentially in time. Broadly speaking, the objective of this paper is to devise polynomial-time algorithms 
for the SCP with probabilistic optimality guarantees, and derive stability conditions for its dynamic counterpart 
(where pickup/delivery requests are generated by an exogenous Poisson process and that serves as a model for DRT 
systems). 

Literature overview. The SCP, being a generalization of the Traveling Salesman Problem, is NP-Hard Q. The 
problem is NP-Hard even on trees, since the Steiner Tree Problem can be reduced to it ||5). In Q, the authors 
present several approximation algorithms for tree graphs with a worst-case performance ratio ranging from 1.5 to 
around 1.21. The 1.5 worst-case algorithm, based on a Steiner tree approximation, runs in linear time. Recently, 
one of the polynomial-time algorithms presented in Q has been shown to provide an optimal solution on almost 
all inputs (with a 4/3-approximation in the worst case) ||6). Even though the problem is NP-hard on general trees, 
the problem is in P on paths [7|. For general graphs, the best approximation ratio is 9/5 and is achieved by an 
algorithm in (8). Finally, an average case analysis of the SCP on trees has been examined for the special case of 
caterpillars as underlying graphs |j9|. 

Dynamic SCPs are generally referred to in the literature as dynamic PDPs with unit-capacity vehicles (1-DPDPs); 
in the dynamic setting pickup/delivery requests are generated by an exogenous Poisson process and the objective 
is to minimize the waiting times of the requests. 1-DPDPs represent effective models for one-way vehicle sharing 
systems, which constitute a leading paradigm for future urban mobility ||3|. They are generally treated as a sequence 
of static subproblems and their performance properties, such as stability conditions, are, in general, not characterized 
analytically. Thorough surveys on heuristics, metaheuristics and online algorithms for 1-DPDPs can be found in 
Pol and pT) . Even though these algorithms are quite effective in addressing 1-DPDPs, alone they do not give 
any indication of fundamental limits of performance. To the best of our knowledge, the only analytical studies 



for 1-DPDPs are |12| and p3|. Specifically, in p2) the authors study the single vehicle case of the problem 



under the constraint that pickup and delivery distributions are uniform; in |13) the authors derive bounds for the 
more general case of multiple vehicles and general distributions, however under the quite unrealistic assumption of 
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three-dimensional workspaces and identical distributions of pickup and delivery sites. 

Contributions. In this paper, we embed the SCP within a probability framework where origin/destination pairs are 
identically and independently distributed random variables within an Euclidean environment. Our random model is 
general in the sense that we consider potentially non-uniform distributions of points, with an emphasis on the case 
that the distribution of pickup sites is distinct from that of delivery sites; furthermore, the graph induced by the 
origin/destination pairs does not have any specific restrictions. We refer to this version of the SCP as the stochastic 
SCP. 



The contribution of this paper is threefold. First, we present the SPLICE algorithm, a polynomial-time algorithm 
for the stochastic SCP, which is asymptotically optimal almost surely; that is, except on a set of instances of 
zero probability, the cost of the tour produced by this algorithm approaches the optimal cost as the number n 



of origin/destination pairs goes to infinity. In practice, convergence is very rapid and SPLICE computes solutions 
for the SCP within 5% of the optimal cost for a number of pickup/delivery pairs as low as 100. The SPLICE 
algorithm has complexity of the order 0{n^^'^) (for arbitrarily small positive constant e), where n is the number of 
pickup/delivery pairs. From a technical standpoint, these results leverage a novel connection between the Euclidean 
Bipartite Matching Problem and the theory of random permutations. 

Second, we provide bounds for the cost of the optimal SCP solution, and also for the solution delivered by 



SPLICE (these bounds are asymptotic and hold almost surely). 

Finally, by leveraging the previous results, we derive a necessary and sufficient stability condition for 1-DPDPs 
for the general case of multiple vehicles and possibly different distributions for pickup and delivery sites. We 
show that when such distributions are different, our stability condition presents an additional (somewhat surprising) 
term compared to stability conditions for traditional spatially-distributed queueing systems. This stability condition 
depends only on the workspace geometry, the stochastic distributions of pickup and delivery points, the demand 
arrival rate, and the number of vehicles. 

Structure of the paper. This paper is structured as follows. In Section [ll] we provide some background on the 
Euclidean Bipartite Matching Problem and on some notions in probability theory and transportation theory. In 
Section III we rigorously state the stochastic SCP, the 1-DPDP, and the objectives of the paper; in Section IV we 



introduce and analyze the SPLICE algorithm, a polynomial-time, asymptotically optimal algorithm for the stochastic 
SCP. In Section [V] we derive a set of analytical bounds on the cost of a stochastic SCP tour, and in Section [Vl| we 
use our results to obtain a general necessary and sufficient condition for the existence of stable routing policies for 



1-DPDPs. Then, in Section VII we present simulation results corroborating our findings. Finally, in Section VIII 



we draw some conclusions and discuss some directions for future work. 

II. Background Material 

In this section we summarize the background material used in the paper. Specifically, we review some results in 
permutation theory, the stochastic Euclidean Bipartite Matching Problem (EBMP), a related concept in transportation 
theory, and a generaUzed Law of Large Numbers. 
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(1) 



A. Permutations 

A permutation is a rearrangement of the elements of an ordered set S according to a bijective correspondence 
cr : 5 5. As an example, a particular permutation of the set (1, 2,3,4) is cr(l) = 3, a{2) — 1, (t(3) = 2, and 
cr(4) = 4, which leads to the reordered set (3, 1, 2, 4). The number of distinct permutations on a set of n elements 
is given by n\ (factorial). We denote the set of permutations over the ji-element ordered set (1, . . . ,7i) by n„. A 
permutation can be conveniently represented in a two-line notation, where one lists the elements of S in the first 
row and their images in the second row, with the property that a first-row element and its image are in the same 
column. For the previous example, one would write: 

12 3 4 
3 12 4 

The identity permutation maps every element of a set S to itself and will be denoted by ai. We will use the following 
elementary properties of permutations, which follow from the fact that permutations are bijective correspondences: 
Prop. 1: Given two permutations cr, o- e n„, the composition aa is again a permutation. 

Prop. 2: Each permutation a £ n„ has an inverse permutation a^^, with the property that a{x) — y if and only 

if cr^^(y) = X. (Thus, note that a^^a = <ti.) 
Prop. 3: For any & E n„, it holds n„ = {o'er, a G n„}; in other words, for a given permutation a playing the 

role of basis, n„ can be expressed in terms of composed permutations. 
A permutation a e n„ is said to have a cycle £ C 5 if the objects in L form an orbit under the sequence 



a{lt)^lt+i for < = 1, . . . ,T- 1, anda{lT) = li, 

where k E C for all t and T is the orbit size (a natural number); given a permutation a, the partition of S into 
disjoint cycles is uniquely determined apart from cyclic reordering of the elements within each cycle (see Figure [TJ. 
Henceforth, we denote by N{a) the number of distinct cycles of <j. In the example in equation ([T]), there are two 
cycles, namely (1,3,2), which corresponds to cr(l) = 3, (7(3) = 2, (t(2) = 1, and (4), which corresponds to 
cr(4) = 4 (see Figure [TJ. 

Suppose that all elements of n„ are assigned probability 1/n!, i.e., 

P[o-] := P[One selects ct] ^, for all a e n„. 

Let Nn denote the number of cycles of a random permutation with the above probability assignment. It is shown 
in p4) that the number of cycles iV„ has expectation E [N^] — log(n)+0(l) and variance var(iV„) = log(n)+0(l); 
here log denotes the natural logarithm. 
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12 3 4 



2 4 



Cycle 1 : {1,3,2} 
Cycle 2 : {4} 



Fig. 1. The two cycles corresponding to tlie permutation: cr(l) = 3, cr{2) = 1, <t(3) = 2, and (7(4) = 4. Cycle 1 can equivalently be 
expressed as (2, 1, 3) or (3, 2, 1). Apart from this cyclic reordering, the decomposition into disjoint cycles is unique. 



B. The Euclidean Bipartite Matching Problem 

Let Xn — {a^i, . . . , x„} and 3^„ — {j/i, . . . , ?/„} be two sets of points in W^. The Euclidean Bipartite Matching 
Problem is to find a permutation a* E n„ (not necessarily unique) such that the sum of the Euchdean distances 
between the matched pairs {{yi,Xcr*{i)) for i = 1, . . . ,n} is minimized, i.e.: 

n n 

- Vt\\ = mill V 11x^(4) - yi\\, 

crGli,, 

1=1 1=1 

where |1 • || denotes the Euclidean norm and n„ denotes the set of all permutations over n elements. Let Qn ■= 
{Xn,yn)\ we refer to the left-hand side in the above equation as the optimal bipartite matching cost Lfa{Qn)', we 
refer to luiQu) ■— Lyi{Qn)/n as the average match cost. 

The Euclidean bipartite matching problem is generally solved by the "Hungarian" method p3) , which runs in 
0{'n?) time. The 0(11?) barrier was indeed broken by Agarwal et al. p6) , who presented a class of algorithms run- 
ning in 0(n^+^), where e is an arbitrarily small positive constant. Additionally, there are also several approximation 
algorithms: among others, the algorithm presented in p7) produces a 0(log(l/e)) optimal solution in expected 
runtime 0(n^+'), where, again, e is an arbitrarily small positive constant. 

The EBMP over random sets of points enjoys some remarkable properties. Specifically, let Xn = {^1, . . . , X„} 
be a set of n points in a compact set Vl C W^, > 3, that are independently, identically distributed (i.i.d.) according 
to a probability distribution with density : il — > M>o; let Xi = {^ii • ■ • j ^n} be a set of n points in Vl that are 
i.i.d. according to the same probability distribution. In |T8j it is shown that there exists a constant /3m. d such that 
the optimal bipartite matching cost Ly[{Qn) = mincren,! J27=i \\-^cr{i) ^ has limit behavior 

1™ = /3^., f <^(x)i-i/<i dx, (2) 

almost surely, where ^ is the density of the absolutely continuous part of the point distribution. The constant /?m.3 
has been estimated numerically as /3m, 3 ~ 0.7080 ± 0.0002 p9) . 

In the case d = 2 (i.e., the planar case) the following weaker result pO) holds with high probability as n +00 
(i.e. with probability 1 — o(l)): 

iM(Q„)/(nlogn)i/2 <^ (3) 
for some positive constant 7. If the probability distribution is uniform, it also holds with high probability that 



LuiQn) /{n log n)^^'^ is bounded below by a positive constant |21 1. 
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To the best of our knowledge, there have been no similar results in the literature that apply when the distributions 
of X points is different from the distributions of y points (which is the typical case for transportation systems). 

C. Euclidean Wasserstein distance 

As noted, little is known about the growth order of the EBMP matching when the distributions of X points is 
different from the distributions of y points. One of the main contributions of the paper is to extend the results 



in 1 18 1 to this general case, for which the following notion of transportation complexity will prove useful. 

Let (pi and Lp2 be two probability densities over J7 C M*^. The Euclidean Wasserstein distance between cpi and 
(p2 is defined as 



Wi(pi,ip2)= inf / \\y- x\\ d'y{x,y), (4) 
7er(ipi,ip2) Jx,y£n 

where V{Lpi^Lp2) denotes the set of measures over the product space $7x17 having marginal densities tpi and if2, 
respectively. The Euclidean Wasserstein distance is a continuous version of the so-called Earth Mover's distance; 
properties of the generalized version are discussed in [22J . 

D. The Strong Law of Absolute Differences 

The last bit of background is a slightly more general version of the well-known Strong Law of Large Numbers 
(SLLN). Let Xi, . . . , Xn be a sequence of scalar random variables that are i.i.d. with mean EX and finite variance. 
Then the sequence of cumulative sums Sn — -^i has the property (discussed, e.g., in |j23j) that 

lim — — — 0, almost surely, 

for any a> 1/2. Note that the SLLN is the special case where a = 1. 

III. Problem Statement 

In this section, we first rigorously state the two routing problems that will be the subject of this paper, and then 
we state our objectives. 

A. The Euclidean Stacker Crane Problem 

Let Xn = {xi, . . . , Xn} and 3^„ = {j/i, . . . , y„} be two sets of points in the d-dimensional Euclidean space R'^, 
where d > 2. The Euclidean Stacker Crane Problem (ESCP) is to find a minimum-length tour through the points 
in Xn U yn with the property that each point Xi (which we call the ith pickup) is immediately followed by the 
point yi (the ith delivery); in other words, the pair {xi,yi) must be visited in consecutive order (see Figure |2|. 
The length of a tour is the sum of all Euclidean distances along the tour. We will refer to any feasible tour (i.e., 
one satisfying the pickup-to-delivery constraints) as a stacker crane tour, and to the minimum-length such tour as 
the optimal stacker crane tour. Note that the ESCP is a constrained version of the well-known Euclidean Traveling 
Salesman Problem. 
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(a) Six pickup/delivery pairs are generated in (b) Daslied arrows combined with tlie solid 
tlie Euclidean plane. aiTows represent a stacker crane tour. 

Fig. 2. Example of Euclidean Stacker Crane Problem in two dimensions. Solid and dashed circles denote pickup and delivery points, respectively; 
solid arrows denote the routes from pickup points to their delivery points. 

In this paper we focus on a stochastic version of the ESCP. Let Xn = {^i, • • ■ tX^] be a set of points in a 
compact set 57 C M'' that are i.i.d. according to a distribution with density ipp : Vl ^ IR>o; let = {Yi, . . . , Yn] 
be a set of points in il that are i.i.d. according to a distribution with density (^d : ^ ]R>o. To obtain the 
relevant transportation problem, we interpret each pair [Xi^Yi) as the pickup and delivery sites, respectively, of 
some transportation demand, and we seek to determine the cost of an optimal stacker crane tour through all points. 
We will refer to this stochastic version of ESCP as ESCP(n, </jp, </jd), and we will write 3^„ ^ ESCP(n, (^p, (^d) 
to mean that Xn contains n pickup sites i.i.d. with density fp, and 3^„ contains n delivery sites i.i.d. with density 
ifY^. An important contribution of this paper will be to characterize the behavior of the optimal stacker crane tour 
through Xn and 3^„ as a function of the parameters n, Lpp, and (^d; throughout the paper we will assume that 
densities ip-p and i^d are absolutely continuous. 

Despite a close relation between the stochastic ESCP and the stochastic EBMP, the asymptotic cost of a stochastic 
ESCP has not been characterized to date. 

B. Dynamic Pickup Delivery Problems with Unit-Capacity Vehicles 

The 1-DPDP is defined as follows. A total of m vehicles travel at unit velocity within a workspace Q,; the 
vehicles have unlimited range and unit capacity (i.e., they can transport at most one demand at a time). Demands 
are generated according to a time-invariant Poisson process, with time intensity A e M>o. A newly arrived demand 
has an associated pickup location which is independent and identically distributed in il according to a density (^p. 
Each demand must be transported from its pickup location to its delivery location, at which time it is removed from 
the system. The delivery locations are also i.i.d. in Vl according to a density tp^). A policy for routing the vehicles 
is said to be stabilizing if the expected number of demands in the system remains uniformly bounded at all times; 
the objective is to find a stabilizing and causal routing policy that minimizes the asymptotic expected waiting times 
(i.e., the elapsed time between the arrival of a demand and its delivery) of demands. 
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This problem has been studied in p3] under the restrictive assumptions (^d = <y9p := and d > 3; in that paper, 
it has been shown that if one defines the "load factor" as 

X¥.^\\Y - X\\ /m, 

where Y and X are two independent random points in f2 with a distribution of density (f, then the condition g <\ 
is necessary and sufficient for a stabilizing policy to exist. However, that analysis — and indeed the result itself — is 
no longer valid if (^d 7^ fp- This paper will show how the definition of load factor has to be modified for the more 
realistic case (/3d 7^ "ySp- Pivotal in our approach is to characterize, with almost sure analytical bounds, the scaling 
of the optimal solution of ESCP(n, (^p, (^d) with respect to the problem size. 



C. Objectives of the Paper 

In this paper we aim at solving the following three problems: 

PI Find a polynomial-time algorithm A for the ESCP which is asymptotically optimal almost surely, i.e., 

lim LA{n)/L*{n) = 1, 

where n is the size (number of demands) of the stochastic instance, is the length of the stacker 

crane tour produced by algorithm A, and L*{n) is the length of the optimal stacker crane tour. 
P2 For the general case ip-£, 7^ cpp, characterize the growth (with respect to the problem size) in the cost of 

the ESCP with almost sure analytical bounds. 
P3 Find a necessary and sufficient condition for the existence of stabilizing policies for the 1-DPDP as a 

function of the problem parameters, i.e.. A, m, ipp, (p-^,, ft. 
The solutions to the above three problems collectively lead to a new class of robust, polynomial-time, provably- 
efficient algorithms for vehicle routing in large-scale transportation systems. 

IV. An Asymptotically Optimal Polynomial-Time Algorithm for the Stochastic ESCP 
In this section we present an asymptotically optimal, polynomial-time algorithm for the stochastic ESCP, which 



we call SPLICE The key idea behind SPLICE is to connect the tour from delivery sites back to pickup sites in 



accordance with an optimal bipartite matching between the sets of pickup and delivery sites. Unfortunately, this 



procedure is likely to generate a certain number of disconnected subtours (see Figure 3(b) 1, and so, in general, the 
result is not a stacker crane tour. The key property we prove is that the number of disconnected subtours grows 
quite slowly, i.e., sublinearly, with the size of the problem. Then, by using a greedy algorithm to connect such 
subtours, one obtains an asymptotically optimal solution to the ESCP with a polynomial number of operations 
(since an optimal matching can be computed in polynomial time). 
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A. The \SPLICE\ Algoritht7i 

The algorithm [splice] is described in pseudo-code. In line[T] the algorithm A4 is any algorithm that computes 
optimal bipartite matchings. After the pickup-to-delivery links and the optimal bipartite matching links are added 
(lines 2]|3i, there might be a number of disconnected subtours (they do, however, satisfy the pickup-to-delivery 
contraints). In that case (i.e., when N > 1), links between subtours are added, e.g., by using a nearest-neighbor 
rul^ Figure [s] shows a sample execution of the algorithm; we refer to the delivery-to-pickup links added in lines 
[T2] and [15] (the green links in Figure ]3]l as connecting links, since they connect the subtours. The complexity 
of SPLICE is dominated by the construction of the optimal bipartite matching, which takes time 0(n^+'^). 



Algorithm SPLICE 

Input: a set of demands S = {(2^1, yi), . . . , (a;„, ?;„ )}, n>l. 
Output: a Stacker Crane tour through S. 

1: initialize a ^ solution to Euclidean bipartite matching problem between sets X — {xi, . . . ,Xn} and y = 
{j/i, . . . , y„} computed by using a bipartite matching algorithm Ai. 

2: Add the n pickup-to-delivery links Xi ^ yi, i = 1, . . . ,n. 

3: Add the n matching Unks yi — ^ a;a(i). i = 1, ■ ■ ■ ,n. 

4: N ^ number of disconnected subtours. 

5: if > 1 then 

6: Arbitrarily order the N subtours Sj, j = 1, . . . ,n, into an ordered set S :— {Si, . . . ,Sn}- 
1: base -i— index of an arbitrary delivery site in Si. 
8: prev -s— base. 
9: for fc = 1 ^ AT - 1 do 
10: Remove link j/prev -> 2: cr(prev)- 

11: next ^ index of pickup site in Sk+i that is closest to 2/prev 
12: Add link j/p^ev ^ a;nexf 
13: prev -S— (T^"' (next). 
14: end for 

15: Add link J/p^ev a;<T(base)- 

16: end if 



' In this paper we use a simple nearest-neighbor heuristic (in 



SPLICE 



lines 



s y to [lej 



for adding connecting links to form an SCP tour 



However, the results of this paper do not depend on this choice, and any connecting heuristic can be used. 
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(c) Line |10| Algorithm state: prev = base = (d) Line |15| Algorithm state: prev = 5. base 

3, fc = 1. The link yg — > xi is removed, next = 3. The link j/5 xi is added an the tour is 

is assigned the value 6. the link y:j — > 2:5 is completed, 
added, prev is assigned the value 5. 



Fig. 3. Sample execution of the |SPLICE| algorithm. The solution to the EBMP is ct(1) = 2, ct{2) = 3, cr(3) = 1, a{i) = 5, cr(5) = 6, 
and cr{6) = 4. Demands are labeled with integers. Pickup and delivery sites are represented by solid and dashed circles, respectively. Pickup- 
to-delivery links are shown as black arrows. Matching links are dark dashed arrows. Subtour connections are shown as lighter, dashed arrows. 
The resulting tour isl-s>2-s>3^6-s>4-s>5-!>l. 



B. Asymptotic Optimality of SPLICE 



In general the SPLICE algorithm produces a number of connecting links between disconnected subtours (i.e., in 
general iV > 1; see Figure [3]). Thus a first step in proving asymptotic optimality of the SPLICE algorithm is to 
characterize the growth order for the number of subtours with respect to n, the size of the problem instance, and 
thus the cost of the extra connecting links. First we observe an equivalence between the number of subtours N 
produced by line |3] and the number of cycles for the permutation a in line [T] 

Lemma 4.1 (Permutation cycles and subtours): The number N of subtours produced by the SPLICE| algorithm 
in line [3] is equal to N{a), where N{a) is the number of cycles of the permutation a computed in line[T| 
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Proof: Let V^. be the set of delivery sites for subtour /c (fc = 1, . . . ,iV). By construction, the indices in 
constitute a cycle of the permutation a. For example, in Figure |3] the indices of the delivery sites in the subtour 
2^1 — > yi — > 2:2 — > 2/2 2:3 — > 2/3 — > xi are {1, 2, 3}, and they constitute a cycle for a since (t(1) ~ 2, (t(2) = 3, 
and (t(3) = 1. Since the subtours are disconnected, and every index is contained by some subtour, then the sets 
Vk (k — 1, . . . , N) represent a partition of {1, . . . , n} into the disjoint cycles of a. This implies that the number 
of subtours N is equal to N{(t). ■ 
By the lemma above, characterizing the number of subtours generated during the execution of the algorithm is 
equivalent to characterizing the number of cycles for the permutation <t. Leveraging the i.i.d. structure in our problem 
setup, one can argue intuitively that all permutations should be equiprobable. In fact, the statement withstands 
rigorous proof. 

Lemma 4.2 ( Equip robability of permutations): Let Qn — 3^„) be a random instance of the EBMP, where 
Xn,yn ^ ESCP(n, (fp,ipu). Then 

¥[(t] = ^ for all a eU, 

where F[a] denotes the probability that an optimal bipartite matching algorithm M produces as a result the 
permutation a. 

Proof: See Appendix |A] for the proof. ■ 
Lemmas 4.1 and 4.2 allow us to apply the result in Section II-A| to characterize the growth order for the number 



of subtours; in particular, we observe that N — N{a) — log(n) + 0(1). For the proof of optimality for SPLICE 
we would like to make a slightly stronger statement: 

Lemma 4.3 (Asymptotic number of subtours): Let d > 2, and let Qn — 3^,i) be a random instance of the 



EBMP, where Xn,yn ^ ESCP(n, tpp, ip-o)- Let iV„ be the number of subtours generated by the SPLICE algorithm 
on the problem instance Qn- Then 

lim Nn/n = 0, 

n— f +00 

almost surely. 

Proof: See Appendix [B] for the proof. ■ 
Remark 4.4: For d > 3, one can similarly prove that limn^^^ Nn/n^~^^'^ = almost surely. 



Having characterized the number of subtours generated by SPLICE we are now ready to prove the main result 
of the section, i.e., the asymptotic optimality of the algorithm. 

Theorem 4.5: Let d > 2. Let Xn be a set of points {Xi, . . . ,Xn} that are i.i.d. in a compact set C M'* and 
distributed according to a density ipp; let 3^„ be a set of points {Yi, . . . , Yn} that are i.i.d. in a compact set J7 C M'' 
and distributed according to a density (po- Then 

lim ^r^'^^f ^ — Ij almost surely. 

n-)-+oo L*[n) 
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Proof: Let Q„ = 3^„). The length of the optimal stacker crane tour through the pickup points and the 
delivery points is bounded below by 

n 

L*{n)>Y,\%~X,\\^Lu{Qn)- (5) 



On the other hand, the number of connecting links added by the SPLICE algorithm is bounded above by the 
number of subtours Nn of the optimal bipartite matching, and the length of any connecting link is bounded above 

by max^ ll^; — ?;||. Hence, ispucEl"-) can be bounded above by 

n 

isPLiCE(«) < \\Yi - X^W + LuiQn) + max ||a; - y\\ Nn 

i—1 

< L* (n) + max \\x — y\\ Nn- 

By the Strong Law of Large Numbers, lim„^-|_oo J2i 11^* ~ ^dl ^ H^ippvdII^ ^ ^11 almost surely. Hence, L*{n) 
has linear growth. Since lim„_j.+oo Nn/n = (by Lemma [43]l, one obtains the claim. ■ 



V. Analytical Bounds on the Cost of the ESCP 

In this section we derive analytical bounds on the cost of the optimal stacker crane tour The resulting bounds 
are useful for two reasons: (i) they give further insight into the ESCP (and the EBMP), and (ii) they will allow us 
to find a necessary and sufficient stability condition for our model of DRT systems (i.e., for the 1-DPDP). 

The development of these bounds follows from an analysis of the growth order, with respect to the instance size 
n, of the EBMP matching on Qn — 3^n), where Xn,yn ^ ESCP {n, ipp, ip-£,). The main technical challenge 
is to extend the results in p8|, about the length of the matching to the case where ipp and (po are not identical. 



We first derive in Section V-A a lower bound on the length of the EBMP matching for the case ipp ^ ifpy (and 



resulting lower bound for the ESCP); then in Section V-B we find the corresponding upper bounds. 



A. A Lower Bound on the Length of the ESCP 

In the rest of the paper, we let C = {C^, . . . , Cl'''} denote some finite partition of Euclidean environment f2 into 
\C\ cells. We denote by 

(/3p(C*) :— I ipp{x)dx 

the measure of cell C* under the pickup distribution (with density (pp), i.e., the probability that a particular pickup 
X is in the ith cell. Similarly, we denote by 



<<5d(C') := / ipD{y)dy 

the cell's measure under the delivery distribution (with density ipp,), i.e., the probability that a particular delivery Y 
is in the ith cell. Most of the results of the paper are valid for arbitrary partitions of the environment; however, for 
some of the more delicate analysis we will refer to the following particular construction. Without loss of generality, 
we assume that the environment C M'' is a hyper-cube with side-length L. For some integer r > 1, we construct 
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a partition Cr of il by slicing the hyper-cube into a grid of r"^ smaller cubes, each length L/r on a side; inclusion 
of subscript r in our notation will make the construction explicit. The ordering of cells in Cr is arbitrary. 

Our first result bounds the average match length luiQn) asymptotically from below. In preparation for this result 
we present Problem [T| a linear optimization problem whose solution maps partitions to real numbers. 

Problem 1 (Optimistic "rebalancing"): 



Minimize >^ 



aij mm \\x - y\\ 



subject to ^ttij = (y9D(C") for all G C, 

3 

a„- = </Jp (C^ ) for all & G C. 

i 

We denote by T(C) the feasible set of Problemjlj and we refer to a feasible solution A{C) := [atj] as a transportation 
matrix. We denote by A{C) :— [a^j] the optimal solution of Problem |Tj which we refer to as the optimistic matrix 
of partition C, and we denote by 1{C) the cost of the optimal solution. 

Lemma 5.1 (Lower bound on the cost of EBMP): Let A'„,3^„ ^ ESCP(7i, (pp, (po), and let Qn ~ 3^„). For 
any finite partition C of J7, liminf„^oo ^M(Qn) > almost surely, where ly^{Qn) is the average length of a 
match in the optimal bipartite matching, and 1(C) denotes the value of Problem [T| 

Proof: See Appendix |B] for the proof. ■ 

We are interested in the tightest possible lower bound, and so we define / := sup^ i(C). Remarkably, the supremum 
lower bound I is equivalent to the Wasserstein distance between c^d and (^p, and so we can refine Lemma 5.1 as 
follows. 

Lemma 5.2 (Best lower bound on the cost of EBMP): Let A",;, ^ ESCP(n, (/3p, ^pn), and let Qn — {Xn, 3^n). 
Then 

lim iid Im{Q n) ^ W((pD,(pp), almost surely. (6) 

n— f oo 

Proof: The lemma is proved by showing that sup^ 1{C) — W((pD, fp)- See Appendix [b| ■ 
Henceforth in the paper, we will abandon the notation I in favor of W((p-D, fp) to denote this lower bound. This 
connection to the Wasserstein distance yields the following perhaps surprising result. 

Proposition 5.3: The supremum lower bound W{(p]:,, fp) of Q is equal to zero if and only if ip-^, ~ Lpp. 
Proof: The proposition follows immediately from the fact that the Wasserstein distance is known to satisfy the 
axioms of a metric on T((fpi,Lpp). Nevertheless, we provide a short alternative proof. 

The proof of the forward direction is by construction: Suppose fi = ^2 — f, let 7 be the measure such 
that 7(J) = Ix:ix,x)ej'fi^^'>'^^ any J G f7 x f7; clearly 7 G T{ip,ip). Then J^y^^\\y~x\\d'y{x,y) = 

The proof of the reverse direction is by contradiction: Suppose (pi ^ Lp2- Then one can choose e > sufficiently 
small, and regions Ai and A2 (where Ai C A2 ^ fl), so that (pi{Ai) > <y92(^2) and — y\\ > e for all x G Ai 
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and y ^ Then for any 7 £ ^{ifii, ^P2) 

\\x - y\\ d-fix, y) > ej {{{x, y) : x e Ai,y ^ A2}) 

= e[(pi(Ai) -7({(a:,y) : x e A^.y e A2})] 
> e[(pi(Ai) -(^2(^2)] > 0. 

■ 

The intuition behind the alternative proof is that if some fixed area A in the environment has unequal proportions 
of X points versus y points, then a positive fraction of the matches associated with A (a positive fraction of all 
matches) must have endpoints outside of A, i.e., at positive distance. Such an area can be identified whenever 



Thus the implication of Lemma 5.2 is that the average match length is asymptotically no less than some constant 
which depends only on the workspace geometry and the spatial distribution of pickup and delivery points; moreover, 
that constant is generally non-zero. We are now in a position to state the main result of this section. 

Theorem 5.4 (Lower bound on the cost of ESCP): Let L*{n) be the length of the optimal stacker crane tour 
through Xn, yn ^ ESCP(n, ipp, (p-o), for compact e M'^, where d > 2. Then 

liminf L*(n)/n > E^p^„ ||r - X\\ + W{vd, <^p), (7) 

almost surely. 

Proof: A stacker crane tour is composed of pickup-to-deUvery links and delivery-to-pickup links. The latter 
describe some bipartite matching having cost no less than the optimal cost for the EBMP. Thus, one can write 

n 

L*{n)/n> - V||y,-X,|| + -LM(Q„). 
n ^ — ' n 

i=l 

By the Strong Law of Large Numbers, the first term of the last expression goes to ¥.,^-p^-^ \\Y — X\\ almost surely. 



By Lemma 5.2 the second term is bounded below asymptotically, almost surely, by iy((pD,V'p). ■ 
Remark 5.5 (Lower bound on the cost of the mESCP): The multi-vehicle ESCP (mESCP) consists in finding 
a set of m stacker crane tours such that all pickup-delivery pairs are visited exactly once and the total cost is 
minimized. The mECSP arises when more than one vehicle is available for service. It is straightforward to show 
that the lower bound in Theorem |5.4| is also a valid lower bound for the optimal cost of the mESCP, for any rn. 



B. An Upper Bound on the Length of the ESCP 

In this section we produce a sequence that bounds Ly[{Qn) asymptotically from above, and matches the linear 
scaling of (|7]). The bound relies on the performance of Algorithm |2j a randomized algorithm for the stochastic 
EBMP. The idea of Algorithm [2] is that each point y & y randomly generates an associated shadow site X', so 
that the collection X' of shadow sites "looks like" the set of actual pickup sites. An optimal matching is produced 
between X' and X which assists in the matching between y and X; specifically, if a; G A" is the point matched to 
X\ then the matching produced by Algorithm |2] contains {y,x). An illustrative diagram can be found in Figure |4] 
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Algorithm |2] is specifically designed to have two important properties for random sets (5„: First, E — y|| is 
predictably controlled by "tuning" inputs — a partition C of the environment and "policy matrix" A{C) — chosen as 
a function of n; second, Lyi{{X' , X)) / n 0+ as n — ^ +00. Later we will show that C and A{C) can be chosen 
so that E \\X' — Y\\ — > W{ipY),ip'p) (as n +00) leading to a bipaitite matching algorithm whose performance 
matches the lower bound of (|2|. 

Algorithm 2 Randomized EBMP (parameterized) 

Input: pickup points X = {xi, . . . , a;„}, delivery points y = {yi, . . . , y„}, probability densities (pp{-) and f-oi'), 

partition C of the workspace, and matrix A{C) E T{C). 
Output: a bi-partite matching between y and X. 
1: initialize X' ^ 0. 

2: initialize matchings M ^ 0; M ^ 0; M ^ 0. 

3: // generate "shadow pickups" 
4: tor y ey do 

5: Let be the cell containing y. 

6: Sample j; j ~ j' with probabihty aij> /ipY){C''). 

7: Sample X' with pdf ipp{ ■ \X' e C^). 

8: Insert X' into X' and {y,X') into M. 

9: end for 

10: M ^ an optimal EBMP between A" and X. 
11: // construct the matching 
12: for X' € X' do 

13: Let {y,X') and be the matches in M and M, respectively, whose ^Y'-endpoints are X' . 

14: Insert {y, x) into M. 
15: end for 
16: return M 



We present the first two properties as formal lemmas: 

Lemma 5.6 (Similarity of X' to X): Let Xi, . . . , Xn be a set of points that are i.i.d. with density ipp; let 
Yi, . . . , y„ be a set of points that are i.i.d. with density (pD. Then Algorithm |2] generates shadow sites X'^, . . . , X'^, 
which are (i) jointly independent of Xi, . . . , X^, and (ii) mutually i.i.d., with density (pp. 



Proof: Lemma 5.6 relies on basic laws of probability, and its proof is relegated to Appendix |B| 



The importance of this lemma is that it allows us to apply equation (|2]l of Section II-B to characterize Lm((X' , X)). 



Lemma 5.7 (Delivery-to-Pickup Lengths): Let F be a random point with probability density cpp,; let X' be the 
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Fig. 4. Algorithm |2] Demands are labeled with integers. Pickup and delivery sites are represented by solid and dashed circles, respectively. 
Pickup-to-delivery links are shown as black arrows. Shadow pickups are shown as dashed squares, with undirected links to their generators 
(delivery sites); also shown are optimal matching links between shadows and pickups. Dashed aiTows show the resulting induced matching. 
Note, this solution produces two discoimected subtours (1,2,3) and (4). 



shadow site of y generated by lines |5]|7] of Algorithm |2] running with inputs C and A{C). Then 

E||X'-y|| < Va„- max \\x - y\\ . 

Proof: Again, see Appendix [B] ■ 

Given a finite partition C, it should be desirable to choose A{C) in order to optimize the performance of 
Algorithm [2j that is, minimize the expected length of the matching produced. We can minimize at least the bound 
using the solution of Problem |2] shown below. We refer to its solution A{C) as the pessimistic matrix 



of Lemma 



5.7 



of partition C. 

Problem 2 (Pessimistic "rebalancing"): 



Minimize a^j max — yll 



subject to ^ttij = (/jd(C*) for all (7* S C, 

i 

a„- = ifip [0) for all £ C. 

i 

Now we present Algorithm [s] described in pseudo-code, which compute^ a specific partition C, and then invokes 
Algorithm [2] with inputs C and A(C). 

Lemma 5.8 (Granularity of Algorithm^: Let r be the resolution parameter, and Cr the resulting grid-based 
partition, used by Algorithm [s] Let y be a random variable with probability density ip-o, and let X' be the shadow 
site of y ~Y generated by lines |5]|7]of Algorithm|2] running under Algorithm|3] Then E \\X' ~ Y\\-W{ipD, (fip) < 
2LVd/r. 

^In the definition of the algorithm, we use the "small omega" notation, where /(■) £ uj{g{-)) implies lim„_>oo f(n)/g{n) = 00. 
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Algorithm 3 Randomized EBMP 

Input: pickup points X — {xi, . . . , Xn}, delivery points y = {yi, . . . , and prob. densities (pp{-) and (pu{-)- 
Output: a bi-partite matching between y and X. 

Require: an arbitrary resolution function res(n) e (jj{n^^'^), where d is the dimension of the space. 
1: r ^ res(n). 

2: C <~ grid partition Cr, of r'^ cubes. 
3: ^ -s- A{C), the solution of Problem [2] 

4: Run Algorithm [2] on {X,y,ipp,ipY),C,A), producing matching M. 
5: return M 



Proof: See Appendix |B] for the proof. ■ 
We are now in a position to present an upper bound on the cost of the optimal EBMP matching that holds in 
the general case when ipp tpo- 

Lemma 5.9 (Upper bound on the cost of EBMP): Let A'„,3^„ ~ ESCP(n, v?p, v?d)j and let Qn = 3^n). For 
d > 3, 

rO., — n.W( in-r\ . in-D^ 

(8) 



hmsup -TTJd - ^y'-PP^'-P-D), 



n— f +00 



almost surely, where 



For d = 2, 



i^{'^v,'^T>) ■= min 

0e{(^P,ipD} 



/3M,d / dx 
Jn 

LuiQn) - nW{ip-D,ipp) 



< 



1, 



(9) 



(10) 



\/nAogri 

with high probability as n +00, for some positive constant 7. 

Proof: See Appendix |B] for the proof. ■ 
We can leverage this result to derive the main result of this section, which is an asymptotic upper bound for the 
optimal cost of the ESCP. In addition to having the same hnear scaling as our lower bound, the bound also includes 
"next-order" terms. 

Theorem 5.10 (Upper bound on the cost of ESCP): Let Xn,yn ^ ESCP(n, (pp, (^d) be a random instance of 
the ESCP, for compact fi e M'^, where d > 2. Let L* (n) be the length of the optimal stacker crane tour through 
XnUyn- Then, for d > 3, 

L*{n) - n[E^p^„ ||r - X\\ + W^((^d, <^p) 



lim sup ■ 

n— f +C30 



< k(.^p,^d), 



almost surely. For d = 2, 



L*{n)-n E^p^„ \\Y - X\\ +W{ipu,^p) 



\Jn logn 

with high probability as n +cx), for some positive constant 7. 



< 



(11) 



(12) 
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Proof: We first consider the case d > 3. Let isPLiCE('^) be the length of the SCP tour through 3^^ generated 



by SPLICE Let Q„ = (A'„,3^„). One can write 



^splice(^) < 11^2 - Xi\\+Lyi{Qn) + max ||x - y\\ Nn 

i—1 



^ - X,\\ - 7iE^p^„ \\Y-X\\]+ (LuiQn) - nW{ifu,^p) 



max \\x — y\\ 



The following results hold almost surely: The first term of the last expression is o{n^^^^'^) (absolute differences); by 



Lemma 



5.9 



the second term is K{ipp, ipj:,)n^^^^'''+o{n^^^^'^); finally, by Remark 4.4 one has lim„_j.+oo Nn/n^^^^'^ — 
0. Collecting these results, dividing on both sides by n^^^/'', and noting that by definition L*{n) < LspucEin), 
one obtains the claim. The proof for the case d = 2 is almost identical and is omitted. 



VI. Stability Condition for DRT Systems 

In the previous section we presented new asymptotic results for the length of the stochastic EBMP and ESCP. 
We showed convergence to Unearity in the size n of the instance, and characterized next-order growth as well 
(equation ( [TT| and equation ([T2]l). Here we use such new results to derive a necessary and sufficient condition for 
the stability of DRT systems, modeled as DPDPs. 

Let us define the load factor as 

g := A [E^p^„ ||r - X\\ + W{^^, ipp)]/m. (13) 



Note that when (pp, ~ ipp, one has W((y9D,<y5p) = (by Proposition 5.3 i, and the above definition reduces to the 
definition of load factor given in |13| (valid for d > 3 and (^d = <^p)- The following theorem gives a necessary 
and sufficient condition for a stabilizing routing policy to exist. 



Theorem 6.1 (Stability condition for DRT systems): Consider the DPDP defined in Section III which serves as 
a model of DRT systems. Then, the condition < 1 is necessary and sufficient for the existence of stabilizing 
policies. 

Proof of Theorem \6.1\ — Part I: Necessity: 
Consider any causal, stable routing poUcy (since the policy is stable and the arrival process is Poisson, the system 
has renewals and the inter-renewal intervals are finite with probabiUty one). Let A{t) be the number of demand 
arrivals from time (when the first arrival occurs) to time t. Let R{t) be the number of demands in the process of 
receiving service at time t (a demand is in the process of receiving service if a vehicle is traveling toward its pickup 
location or a vehicle is transporting such demand to its delivery location). Finally, let Si be the servicing time of 
the jth demand (this is the time spent by a vehicle to travel to the demand's pickup location and to transport such 
demand to its delivery location). 
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The time average number of demands in the process of receiving service is given by 

9 



1 /■* 

lim - / RMdr. 

t->+oo t J^^o 



By following the arguments in |24 page 81-85, Little's Theorem], p can be written as: 



t^+oo t t^+oo A{t) t-^+oo t 

where the first equality holds almost surely and all limits exist almost surely. The second limit on the right is, by 
definition, the arrival rate A. The first limit on the right can be lower bounded as follows: 

t-J-+oo A{t) t-i.+oo A{t) ' 

where L*{A{t)) is the optimal length of the multi-vehicle stacker crane tour through the A{t) demands (i.e., is 



the optimal solution to the multi-vehicle ESCP - see Remark 5.5 ). For any sample function, L* {A{t)) / A{t) runs 



through the same sequence of values with increasing t as L*{n)/n runs through with increasing n. Hence, by 



Theorem 5.4 and Remark 5.5 we can write, almost surely. 



Collecting the above results, we obtain: 

g>X[E^,^^\\Y-X\\+W{ipu,ipp)], 

almost surely. Since the poUcy is stable and the arrivals are Poisson, the per-vehicle time average number of demands 
in the process of receiving service must be strictly less than one, i.e., g/m < 1; this implies that for any causal, 
stable routing policy 

A [E^p^, \\Y - X\\ + W{^j,,^p)]/m < 1, 

and necessity is proven. 



Proof of Theorem 6.1 —Part II: Sufficiency : The proof of sufficiency is constructive in the sense that we design 



a particular policy that is stabilizing. In particular, a gated policy is stabilizing which performs the following steps 



any time all servers are idle: (1) applies algorithm SPLICE to determine tours through the outstanding demands 



(2) splits the tour into m equal length fragments, and (3) assigns a fragment to each vehicle. 

Consider, first, the case d — 2. Similarly as in the proof of Theorem 4.2 in | |25] , we derive a recursive relation 
bounding the expected number of demands in the system at the times when new tours are computed. Specifically, 
let ti, i > 0, be the time instant at which the ith SPLICE tour is constructed (i.e. the previous servicing round was 
completed); we will call this instant epoch i. We refer to the time interval between epoch i and epoch i + 1 as the 
ith iteration; let Ui be the number of demands serviced during the ith iteration (i.e. all outstanding demands at its 
epoch), and let Ci be the interval duration. 
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Demands arrive according to a Poisson process with rate A, so we have E [jii+i] = A E [Ci] for all epochs. The 
interval duration Ci is equal to the time required to service the demands of the ith group of demands. One can 
easily bound E [Ci] as 

E [C,] < E [LspucEin^)]/m + + D(f7), (14) 

where LsPLiCE('^i) denotes the length of the SPLICE tour through the ith iteration demands, and D(il) = max{||p— 
q\\ \ p,q E il} is the diameter of ft; the constant terms account conservatively for (i) extra fragment length incurred 
by splitting the tour between vehicles, and (ii) extra travel required to reach the current tour fragment from the 
endpoint of the previous fragment. Let 

ip{n):=n E^^^^ \\Y - X\\ + W{ipB, ifip) + ^\/n\ogn + o (^^ynlognj , 

and let q{n) be the probability that given n demands equation ([12]) does not hold, i.e.; 



lin) :=!-: 



L 



SPLICE 



(n) < ^{n) 



Now, let £ be an arbitrarily small positive constant. By Theorem 5.10 lim 



q{n) = 0; hence, there exists 



a number h such that for all n > n one has q{n) < e. Then, the length of the optimal stacker crane tour 
through the Ui demands can be upper bounded as follows (where in the third step we use the trivial upper bound 
Lsplice('t-) < 2nD(ri), which is always valid): 



+00 



E [isPLICE('^i)] = X! [-^SPLICE 



Ui — n\ r n,- = n 



n=0 

+00 

E 



E [isPLiCE('^») I (12} holds, Ui = n] P[([T2]) holds] Ui = n] + 



= l-q{n) 



E [LspucEin-i) I ( fT2] l does not hold, Ui = n] ¥\(l2\ does not hold| rii = 



=q(n) 



+00 



+00 

E 



E [ip{ni) I ([12]) holds, Ui = n] (1 - q{n)) + 

E [2 Hi D(fi) I ([T2| does not hold, Ui = n] q{n) 



?A(n)(l-g(n)) + (2nD(f7))(j(n) 

< E [^(n,)] + ^ (2 n q{n) P[n, = 

= E [^(n,)] + XI (2 nT>{9)) q{n) f>[n, = n] + ^ (2 n D(f^)) g( 

n=0 

< E [^l}{n,)\ + D(f7) + e 2 V){Vt) E [n,\. 



+ 00 



<£ 
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Then, one can write the following recurrence relation: 

E [Q] < E [LsPLiCE(n»)]/TO + 2D(r!) 

< E [■il;{n,)]/m, + e 2 E [n,]/m + D{n){n'^/m + 2) 

< E [n,] [(E^p^o ||r - X\\ + W{ipB,^p)) + s 2 /m 

+ E [7 ^naog + o ( log n^)] /to + D(17)(nVm + 2). 
Now, let (5 > 0; theij^ for all x > 1 

X log X + o{\/ X log x) < c{d) + 6 X, 

where c{6) E M>o is a constant. Hence, we can write the following recursive equation for trajectories i i-> E [rii] S 
M>o: 

E[n,+i] = AE [Ci] <E[n,] [g + e \ 2D{fl)/m + \ S /m] + X-f c{5)/m + D{fl){n^ /m + 2). 

Since g < 1 and e and 5 are arbitrarily small constants, one can always find values for e and S, say e and (5, such 
that 

a{g) + e A2D(fi)/TO + A75/TO < 1. 
Hence, for trajectories i i-> E [ui] e M>o we can write, for all i > 0, a recursive upper bound: 

E[n,+i] < a{g)E[ni] + Xjc{S)/m + D{n){n^/m + 2), (15) 

with a{g) < 1, for any g E [0, 1). 

We want to prove that trajectories i 1— > E [rii] are bounded, by studying the recursive upper bound in equation 
([TSj. To this purpose we define an auxiliary system, System-X, whose trajectories i 1-^ Xi E IR>o obey the dynamics: 

Xi+i = a{g) Xi + \-fc{S)/m + 2D(f7)(nV"^ + 1), (16) 

with Xf) = uq. By construction, trajectories i h- > Xi upper bound trajectories i 1-^ E [rii]. One can easily note that 
trajectories Xi are indeed bounded for all initial conditions, since the eigenvalue of System-X, a{g), is strictly 
less than one. Hence, trajectories i i-^ E [rii] are bounded as well and this concludes the proof for case d — 2. 

Case d > 3 is virtually identical to case d = 2, with the only exception that equation ( [T2] i should be replaced 
with equation ( fTTj ), and the sublinear part is given by x^^^/'^ (the fact that for > 3 the inequalities hold almost 
surely does not affect the reasoning behind the proof, since almost sure convergence implies convergence with high 
probability). ■ 



Note that the stability condition in Theorem 6.1 depends on the workspace geometry, the stochastic distributions 



of pickup and delivery points, the demands' arrival rate, and the number of vehicles, and makes explicit the roles 

'Consider, first, tlie case without tlie o(y/x \ogx) term. In tliis case, let c{S) = ^ (^og ^ — 1^- Then, by using Young's inequality, one 
can write: i (log jj — 1) + S x — \/x log x > i (log -p^— l)+5a; — 4p — '^r^ =• '/'(^)- Then, one can easily show that ^(^x) > 



for all z > 1. The case with the o{\/x logx) term is similar and is omitted. 
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of the different parameters in affecting the performance of the overall system. We believe that this characterization 
would be instrumental for a system designer of DRT systems to build business and strategic planning models 
regarding, e.g., fleet sizing. 

Remark 6.2 (Load factor with non-unit velocity): In our model of DPDPs we have assumed, for simplicity, that 



vehicles travel at unit velocity. Indeed, Theorem 6.1 holds also in the general case of vehicles with non-unit velocity 
V, with the only modification that the load factor is now given by 



Q 



A[E^p^j|y-x|| 

V m 



VII. Simulation Results 

In this section, we present simulation results to support each of the theoretical findings of the paper. Our 
simulation experiments examine all three contributions of the paper; in particular, we discuss (i) performance 



of file SPLICE algorithm, (ii) scaling of the length of the EBMP, and (iii) stabilizability of the DPDP. 



A. Performance of SPLICE 



We begin the simulation section with a discussion of the performance of the SPLICE algorithm. We examine 



(i) the rate of convergence of SPLICE to the optimal solution, (ii) the runtime for the SPLICE algorithm, and 



(iii) the ratio between the runtime of SPLICE and that of an exact algorithm. In all simulations we assume that 
the pickup/delivery pairs are generated i.i.d. in a unit cube and that ipp and (pu are both uniform. The Bipartite 
Matching Problem in line [T| of [SPLICE is solved using the GNU Linear Programming Toolkit (GLPK) software on 



a linear program written in MathProg/AMPL; for comparison with SPLICE the Stacker Crane Problem is solved 



exactly using the GLPK software on an integer linear program. (Simulations were run on a laptop computer with 
a 2.66 GHz dual core processor and 2 GB of RAM.) 



Figure 5(a) shows the ratios -Lsplice/^* observed in a variety of randomly generated samples (twenty-five trials 
in each size category). One can see that the ratio is consistently below 20% even for small problem instances 
(n ~ 10) and is reduced to ~ 5% for n > 80. Hence, convergence to the optimal solutions with respect to the 



problem size is fairly rapid. In practice, one could combine SPLICE with an exact algorithm, and let the exact 



algorithm compute the solution if n is less than, say, 50, and let SPLICE compute the solution when n > 50 



Figure 5(b) shows the ratios Tsplice/T* with respect to the size of the problem instance n (the problem instances 
are the same as those in Figure 5(a)| l, where T* is the runtime of an exact algorithm. One can observe that the 
computation of an optimal solution becomes impractical for a number n ~ 100 of origin/destination pairs. 



Finally, Figure 5(c) shows the runtime Tsplice of the SPLICE algorithm with respect to the size of the problem 
instance n (the problem instances are the same as those in figure |5(a)[ ). One can note that even for moderately large 
problem instances (say, n ~ 100) the runtime is below a second. 
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(a) Cost factor for SPLICE as a function of the problem size n. (b) Runtime factor for SPLICE as a function of the problem size 
Each column records 25 random observations. n. Each column records 25 random observations. 
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(c) Runtime of the SPLICE algorithm as a function of the problem 
size n. Each column records 25 random observations. 



Fig. 5. Performance of the SPLICE algorithm over twenty-five trials in each size category. Figure p(a)| shows observed ratios isPLlCfi/^'* 
Figure [5(b)| shows the ratios Tspuce/T* ', Figure [5(c)| shows the runtime Tspdce of the |SPLICE| algorithm. 



B. Euclidean Bipartite Matching — First- and Next- Order Asymptotics 

In this section, we compare the observed scaling of the length of the EBMP as a function of instance size, with 
what is predicted by equations (|7]i and ( fTTj ) of Sections V-A and V-B respectively. We focus our attention on two 
examples of pickup/delivery distributions ((/Sp, (/Sd): 

Case I — Unit Cube Arrangement: In the first case, the pickup site distribution tp-p places one-half of its probability 
uniformly over a unit cube centered along the x-axis at a; = —4, and the other half uniformly over the unit cube 
centered at x = — 2. The delivery site distribution i/jd places one-half of its probability uniformly over the cube at 
X = —4 and the other half over a new unit cube centered at a; = 2. 

Case II — Co-centric Sphere Arrangement: In the second case, pickup sites are uniformly distributed over a sphere 
of radius R = 2, and delivery sites are uniformly distributed over a sphere of radius r = 1. Both spheres are 



February 8, 2012 



DRAFT 



24 



centered at the origin. 



Figures 6(a) and 6(b) show samples of size n = 100, drawn according to the distributions of Case I and Case 
II respectively; the left plots show the samples alone, while the plots on the right include the links of the optimal 
bipartite matching. 

The cases under consideration are examples for which one can compute the constants W (Wasserstein distance) 
and K of equation ( |TT] l exactly. In the interest of brevity, we omit the derivations, and simply present the computed 
values in Table [l] The extra column "^(1^90, (/jp)" of the table shows a new smaller constant that results from 
bringing the min operation inside the integral in equation (|9]l. 







k(vd,'Pp) 




Case I 


2 


0.892 


Si 0.446 


Case II 


0.75 


^ 1.141 


Si 0.285 



TABLE I 

Values computed for the constants VF(ipDi '/'p) and k(ipd, V'p) in equations jTJ and (TT), for Case I and Case II, 

RESPECTIVELY; ALSO K(ipDi V'p) IN EACH CASE, THE RESULT OF BRINGING THE min OPERATION INSIDE THE INTEGRAL IN EQUATION j9). 



The simulation experiment is, for either of the cases above, and for each of seven size categories, to sample 
twenty-five EBMP instances of size n randomly, and compute the optimal matching cost Lm of each. The results of 
the experiment are shown in Figure |7] Figure 7(a) (top) shows a scatter plot of (n, Lyi/n) with one point for each 



trial in Case I; that is, the x-axis denotes the size of the instance, and the j;-axes denotes the average length of a 
match in the optimal matching solution. Additionally, the plot shows a curve (solid line) through the empirical mean 
in each size category, and a dashed line showing the Wasserstein distance between ip-£, and ipp, i.e. the predicted 
asymptotic limit to which the sequence should converge. Figure [7(b)] (top) is analogous to Figure [7(a)| (top), but for 
random samples of Case II. Both plots exhibit the predicted approach of Ly^/n to the constant W{ipY>i'Pp) > 0; 
the convergence in Figure |7(b) (top) appears slower because W is smaller. Figure 7(a) (bottom) shows a scatter 



plot of {n, {Lm ~ W^)/?!^/"^) from the same data, with another solid curve through the empirical mean. Also shown 
are k and k (dashed lines); recall that k is the theoretical asymptotic upper bound for the sequence (equation ([TT]i). 



Figures 7(b) (bottom) is again analogous to Figures 7(a) (bottom), and both plots indicate asymptotic convergence 



to a constant no larger than k(93d,<^p)- In fact, these cases give some credit to a developing conjecture of the 
authors. The conjecture is that the minimization in (|9]) can be moved inside the integral to provide an upper bound 
like equation ([TT]), but with a smaller (often much smaller) constant factor, i.e. K{(p-£, , (pp) . 

C. Stability of the DPDP 

We conclude the simulation section with some heuristic validation of equation ( [T3| l and the resulting threshold 
A* separating stabilizable arrival rates from unstabilizable ones. The main insight of this section is as follows. Let 
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(a) Case I sample (n = 100) with and without optimal matching. 




(b) Case II sample (n = 100) with and without optimal matching. 

Fig. 6. Samples of size n = 100, drawn according to the distributions of Case I (Figure |6(a)) and Case II (Figure |6(b)^ . Pickup sites are 
shown as tiiangle markers; delivery sites are shown as circles. Left plots show the samples alone; right plots include links of the optimal bipartite 
matching. 
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Fig. 7. Scatter plots of (n, Lu/n) (top) and (n, (Lm — W)/n?^^) (bottom), with one point for each of twenty-five trials per size category. 
Figure [7(a)] shows results for random samples under the distribution of Case I; Figure [7(b)] shows results for random samples under the distribution 
of Case II. 



TT be a policy for the DPDP that is perfectly stabilizing, i.e., stabilizing for all A < A*. We consider the system 
(DPDP(A), tt), where A > A*. Clearly, since A > A*, the number of outstanding demands in the system grows 
unbounded. Still, demands arrive at rate A in time average, and we should expect the policy to serve demands at 
an average rate of A* (i.e., the fastest rate under vr). Thus, the number of outstanding demands should grow at an 
average rate of A — A*. Since we can control A in simulation, we can use this insight to estimate A*, e.g., by the 
simple calculation A — n{T)/T after sufficiently large time T, where n{T) is the number of outstanding demands 
at time T. We focus our discussion on the single-vehicle setting, but results for multiple vehicle systems have 
been equally positive. Consider again the cases of Section VII-B Table |ll] shows computed threshold A* for both 
cases — and the statistics essential in computing it — as well as the estimate of A* after time T — 5000. 





E[Y - X] 


W 


A* = {E[Y - x] + wy'^ 


A* estimate after T = 5000 


Case I 


Ki 3.2 


2 


Si 0.190 


0.20 


Case II 


1.66 


0.75 


Si 0.415 


0.42 



TABLE II 

Stabilizability tresholds a* for Case I AND Case II OF Section [ViI-B| with relevant statistics. Also, an estimate of a* in 

EACH CASE AFTER SIMULATION FOR TIME T = 5000. 



Our simulations were of the nearest-neighbor policy (NN); i.e., the vehicle's ith demand is the demand whose 
pickup location was nearest to the vehicle at the time of delivery of the (i — l)th demand. (The simulated rate 
of arrivals A was 1.) Although a proof that the NN policy is perfectly stabilizing is currently not available, it 
has been observed that such policy has good performance for a variety of vehicle routing problems; it also has a 
fast implementation where large numbers of outstanding demands are concerned. In both cases, the estimated and 
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computed A* were quite close (within 5% of each other). 

VIII. Conclusion 

In this paper we have presented the |SPLICE| algorithm, a polynomial-time, asymptotically optimal algorithm for 



the stochastic (Euclidean) SCR We characterized analytically the length of the tour computed by SPLICE and 
we used such characterization to determine a necessary and sufficient condition for the existence of stable routing 
policies for the 1-DPDP, a dynamic version of the stochastic SCP Our results would provide a designer of DRT 
systems with essential information to build business and strategic planning models regarding, e.g., fleet sizing. 

This paper leaves numerous important extensions open for further research. First, we are interested in precisely 
characterizing the convergence rate to the optimal solution, and in addressing the more general case where the pickup 
and delivery locations are statistically correlated. Second, we plan to develop policies for the dynamic version of 
the SCP whose performance is within a constant factor from the optimal one. Third, while in the SCP the servicing 
vehicle is assumed to be omnidirectional (i.e., sharp turns are allowed), we hope to develop approximation algorithms 
for the SCP where the vehicle has differential motion constraints (e.g., bounded curvature), as is typical, for example, 
with unmanned aerial vehicles. In addition to these natural extensions, we hope that the techniques introduced in 
this paper (e.g., coupling the EBMP with the theory of random permutations) may come to bear in other hard 
combinatorial problems. 
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Appendix A 

The Permutation Probability Assignment of Optimal Bipartite Matchings 



In this section we provide the rigorous proof of Lemma 4.2 in Section IV-B showing the equiprobability 
of permutations produced by an optimal bipartite matching algorithm M.. Let Xn — {xi, . . . , a;„} and = 
{yi, . . . ,?/„} be two sets of points in il C M"^, e.g. forming an instance of the EE MP in environment VL. Consider 
s = concat(xi, . . . , a;„, ?/„), a column vector formed by vertical concatenation of xi, yi, . . . , a;„, Note that 
the set ri'^" of such vectors, i.e. the span of the instances of the EBMP, is a full-dimensional subset of M'*'^"'. Let 
n* : fi^" — > 2^^" denote the optimal permutation map that maps a hatch s E $7^" into the set of permutations that 
correspond to optimal bipartite matchings (recall that there might be multiple optimal bipartite matchings). Let us 
denote the set of batches that lead to non-unique optimal bipartite matchings as: 

Z := {se r!2"| |n*(s)| > l}, 

where |n*(s)| is the cardinality of set 11* (s). 

In Lemma [4!2l A4 may be any algorithm that computes an optimal bipartite matching (i.e., a permutation that 
solves an EBMP). According to our definitions, the behavior of such an algorithm can be described as follows: 
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given a batch s e il^" it computes 

{unique a e U* (s) if s e fl^" \ Z, 
some a G n*(s) otherwise. 
Thus, the behavior of a bipartite algorithm on the set Z can vary; on the other hand, we now show that set Z has 
Lebesgue measure zero so that the behavior of an algorithm on this set is immaterial for our analysis. 
Lemma A.l (Measure of multiple solutions): The set Z has Lebesgue measure equal to zero. 

Proof: The strategy of the proof is to show that Z is the subset of a set that has zero Lebesgue measure. 
For cr', a" £ n„, a' ^ a" , let us define the sets: 

n n 



i=l i=l 



let us also define the union of such sets: 



H 1^ Ha', cr"- 

(T',(T"en„ 

The equality constraint in the definition of T-La^a" implies that T-La^a" ^ ]R'^(2n)-i^ which has zero Lebesgue 
measure in M'^^^"). Hence, the Lebesgue measure of Ha', a" is zero. Since % is the union of finitely many sets of 
measure zero, it has zero Lebesgue measure as well. 

We conclude the proof by showing that Z C H. Indeed, if s e Z, there must exist two permutations a' 7^ a" 

such that X;r=i W^a'it) - VtW = miiV E"=i W^-yii) - and X^^Li \\^cr"{i} ~ ViW = miiV ELi \\^<y{i} - VtW^ i-^-, 
there must exist two permutations a' ^ a" such that 

n n 

\\Xa'{i) - ViW = \\^<^"i^) ~ 2^*11' 

i=l 1=1 

which implies that s E H. Hence, Z C H and, therefore, it has zero Lebesgue measure. 



Now we present the proof of Lemma 4.2 which gives the probability that Ai produces as a result the permutation 
cr for X„, ^ ESCP(n, (^p, (^d); we call such probability P[cr]. 

Proof of Lemma \4.2\ We start by observing that it is enough to consider a restricted sample space, namely 



^2n ^ 2. One can write P[s E Z] = J^^^ 'p{s)ds, where Lp{s) denotes the product Y[i=i 'Pp{^i)'PD{yi)- Because 



of our continuity assumptions on probability distributions. Lemma A.l implies P[s E Z] = 0. Thus, by the total 
probability law, 

V[a]=¥[M{s) = a\sEn^"\Z]. (17) 
For each permutation a E n„, let us define the set 

Sa {s E 172"\Z|7W(s) =a} . 
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Collectively, sets 5^ form a partition of fi^" \ Z. This fact, coupled with equation ( [T7] l, implies 

P[cr] =P[s e 5,,]. 

For a permutation a <E n„, let us define the reordering function : ri^" — > 51^" as the function that maps a batch 
s — concat(a;i, . . . , a;„, ?/„) into a batch s' = concat(a;CT(i), 2/1, • • ■ , Xcr(n),yn)- Alternatively, let Ej £ ^dx2nd 
be a block row matrix of 2n dx d square blocks whose elements are equal to zero except the (2j — l)th block that 
is identity; let Fj E ^dx2nd s,ixch a block matrix, but whose elements are all zero except the 2jth block that is 
identity. Then in matrix form the reordering function can be written as ga{s) — Pa s, where is the 2nd x 2nd 
matrix defined by 



Pn ■■= 



Note that | det(Pcr)| = 1 for all permutations ct; also. Prop. 2 of Section II-A implies P^-^ = P^^- We now show 
that ga-{Sa) = for all permutations a £ n„, recalling that ai denotes the identity permutation (i.e., a{i) ~ i 
for i — 1, . . . ,n): 

• gaiSa) C Sa,- Let s e S^- Then by definition Y.'i=i ~ ViW = min(Ten„ 1]"=! \\^<^{i) ~ ViW'- moreover 
(7 is the unique minimizer. We want to show that gais) G Sa^, where gais) has the form 
concat(i;o.(i), yi, . . . ,xs-(n),yn)- Let s — ga{s)', indeed, cri is an optimal matching of s (by inspection), i.e., 
(Ti G n*(s). Suppose, however, there is another optimal matching a ^ such that a E n*(s). Then aa is an 
optimal matching of s (Prop. 1); yet this is a contradiction, because aa ^ a. Therefore, we have that s E Sa^ 
for all s E Sa- 

• C ga{Sa). Let s E S^^. Then by definition Yl7=i W^i ^ y^ll = miiven„ I]"=i \\x<7{i) - yi\\\ moreover 
cTi is the unique minimizer. Note that g^ is an injective function (since the determinant of P^ is nonzero); 
let s be the unique batch such that s ~ ga{s), i.e., s — conca.t{xa-i(i),yi, ■ ■ ■ ,Xa-i(n)Tyn) (Prop. 2). We 
want to show that s E S^- Because \\xi ^ ViW = X]"=i II^ct(ct-MO) ^ o- is an optimal matching of 
s, i.e., a E 11* (s). Suppose there is another optimal matching a ^ a such that a E 11* (s). Again, this is a 
contradiction, since aa~^ ^ a\, and a\ is the unique optimal matching for batch s. We conclude that s E Sa 
for all s E S^^ . 

We are ready to evaluate the probabilities of permutations as follows: For any permutation u we have P[(t] = 
P[s E Sa\ ~ jsi£s- where Lp{s) denotes YTi=i 'Ppixi)f'D{yi)- We use variable substitution s = ga{s) = 

P^s and the property ga{Sa) — S^^, and we apply the rule of integration by substitution: J-^^, ip{s)ds = 
J^^g^ (/?(Pr^s) |det(P^)|~^^ ds. Observing that 

=1 

n 

f{Pa^s) = ip{Pa-is) = Y[fp{xa-i(^))ipD{yi), 

4 = 1 

and that 

n n 

Y[ipp{xa-i{i})(PB{yt) = Y['pp{xi)'fiB{yi) = fis), 
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we obtain 



Combining these results, we conclude P[a] = P[(Ji] for all a E n„, obtaining the lemma. 

Appendix B 
Proofs of Other Lemmas 

A. Lemmas of Section Uv\ 



Proof of Lemma 4.3 For any e > 0, consider the sequence E of events, where 

:y„) : Njn>e} 



Eri 



4.1 



the number of 



or, equivalently, =|(A'„,3^„) : (iV„ - V.Nn) + {V^Nn - log(n)) +log(n) > enj. By Lemma 
disconnected subtours is equal to the number of cycles in the permutation a computed by the matching algorithm M 
in Une[T| Since, by Lemma 4.2 all permutations are equiprobable, the number of cycles has expectation and variance 



both equal to log(n) + 0(1). Therefore, we conclude that Nn has expectation and variance both log(n) + 0(1) 
Hence, we can rewrite the events £"„ as: 



E„ 



--[{Xn.yn) : iV„ - > en + o(n)}. 



Applying Chebyshev's inequality, we obtain (for n' sufficiently large, yet finite) 

n— n—n' L \ J\ 



Since this summation is finite, we can apply the Borel-Cantelli lemma to the sequence of events En and conclude 
that P[limsup„_^_)_oQ £"„] ~ 0. Finally, since e can be chosen arbitrarily small, the upper limit of the claim follows 
(the lower limit holds trivially). ■ 

B. Lemmas of Section [V| 

Proof of Lemma [577} Let a denote the optimal bipartite matching of (5„. For a particular partition C, we define 
random variables :— |{A: : e 0%Xg.(fc) G O-'H /n for every pair {C^,C^) of cells; that is, denotes the 
fraction of matches under a whose 3^-endpoints are in C* and whose ^Y-endpoints are in . Let Tn be the set 
of matrices with entries {a^- > 0}^ ^^-^ such that a^- = n jn for all e C and — 

n (7*1 /n for all C" G C; note {a^} itself is an element of Tn- Then the average match length Zm('3ti) is 
bounded below by 



1 

im^Qu) = - ^ ii^o-(fe) - ^fcii > ' 



_ mm a; — W 

fe=i 11 



> min > a,,,- min llx — ull 
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The key observation is that lini„_j.oo n C^}\ jn — (pp(C^), and lim„^oo |{3^n ^1 C**}] 1"^ = ^'d{C'^\ almost 
surely. Applying standard sensitivity analysis (see Chapter 5 of [26 J), we observe that the final expression converges 
almost surely to 1(C) as n — > +oo. 

Specifically, for any finite partition C and e > 0, consider the sequence of events E, where 



min > a, 



min \\x — y\\ — 1{C) 



> e 



Consider Problem [l| let Ad,; for all i be the dual variables associated with the constraints J2j onj — V'd(C''); let 
Xp.j for all i be the dual variables associated with the constraints aij = fp{C^). (These are all finite constants.) 
Also, let AD,(n) = \yn n C^l /n - (^d(CO for all i; let Apj(n) = \Xn n C^\/n - Lpp{C^) for all j. Through 



sensitivity analysis we obtain 



AeTr 



; mm 



y\\ - 1{C) = J2 Ad^AdzW + o(AD.(n)) + J2 ApjAp,(n) + o(Ap,(n)) 



and so 



mm 

Aetr 



min ||x — y\\ — 1(C) 



< ^IAd^IAd^WI +o(AD.(n)) 

i 

+ ^|Ap,||Ap,(n)|+o(Ap,(n)). 



(18) 



Thus we can choose > for all i and ep^ > for all j so that the right-hand side of ( [TSj l is less than or equal 
to e as long as IAdJ < eoi for all i and |Apj | < epj for all j. Let us now define the events -Bad; for all i and 
-E-Apj for all j, where 

-E^AD.n = {('^«,3^n) : |AD,|>eDj and ^Ap^n = { Xi) : |Apj|>epj}. 

Note that the Strong Law of Large Numbers gives P[limsup„^^^ Ej^^^n] — for all i, and P[lim sup„^_|_o^ -EAp^n] 
for all j. Observing that 

P[limSUp„^+^ En] < ^P[limSUp„^+^ £^An,n] + ^P[limSUp„_^+„^ ^^Ap^-n] = 
i 3 

we obtain the claim. ■ 



Proof of Lemma 5.2- First, we show that I < T4^((pD, V^p)- Fix 7*, some solution arbitrarily close to the 
infimum of Q. Given a partition C let us define the distance function d{xi^X2) = '[ninr^^c(xi),yec{x2) Wu ~ ^W' 
where C( ) maps points to their containing cells; note that d is everywhere a lower bound for the Euclidean distance. 



Noting that the matrix 



ec\x2ec 



J d"f*{xi,X2) 



satisfies the constraints of Problem 1 we have 







W{ipB,'Pp) + e = 



\\X2 - Xi\\ d'J*(xi,X2) 



> 



d{xi,X2) d'J*(^Xi,X2) 



Xi,X2&f^ 



V min lly-xll / d'y*{xi,X2) 



X2<£C^ 



>l(C). 
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Since this inequality holds for all C and e arbitrarily small, we conclude the first part. Next, we show that for any 
e > 0, there exists a partition C so that W{(p-D, (pp) < 1{C) + e. For example, we may choose construction C,. with 
resolution r sufficiently fine so that \\x2 — a;i|| < d{xi,X2) + e everywhere. Solve for A{C), the solution yielding 
1{C). Choose 7 G T{(pY),Lpp) to be any distribution that satisfies /^.^g^i X2eci "^7(^1' ^2) = «y for all i and j (in 
general, there are infinitely many). Then we have 1{C) = J d{xi,X2) dj(xi,X2) > J \\x2 — xi\\ dj{xi,X2) — e > 

M^((/JD, Vp) - e- ■ 

Proof of Lemma 5.6 The delivery sites Yi, . . . , y„ are i.i.d. and independent of Xi, . . . , X„. Thus, under the 
explicitly point-wise independent sampling of Algorithm|2j the n joint variables [Y^, Jk,X'f,) fork— 1, . . . , n, retain 
the same set of independencies. This suffices to prove (i) and the mutual independence of (ii). To complete the proof, 
we derive the marginal distribution of some X' explicitly, by writing the joint distribution of (F, J,X'), and then 
eliminating {Y,J). The joint distribution can be expressed as (pYjx'{y,j,x') = (fuiu) pi'ob(j|2;) fx'\.j{x'\j) = 
Viy{y) X]i=i "yjfcy ^v{x'\X' S C^); here Ic{x) is the indicator function accepting the condition a; e C. 
Eliminating Y, we obtain ipjx'{j,x') = Jy^^tpY.jx'iy,j,x') dy = Jy^^ipB{y) ELi "ypfc-f fpi^'l^' ^ 
CJ) dy = <y9p(a;'|X' e C^) J2f=i = ^pix'\X' e C^)ipp{C^). Finally, taking the sum over j = l,...,r'^, we 
obtain ipx'{x') = (pp{x'), completing the proof. ■ 



Proof of Lemma 5. 7 The proof is by derivation: 

K\\X'~Y\\^ I \\x' ~y\\Y,^YjX'{y.j,x') dx' dy 



x',yeC' 



max lla; — y| 



\\x' -y\\^D{y)vp{x'\X' eC^) dydx' 



L 



^p,{y)ipp{x'\X' £C^) dydx' 



x',yeC' 



: max 

yeCi,xec^ 



\\x - 2/11 



Proof of Lemma 5.8- Algorithm [3] uses partition Cr, and optimal solution A{Cr) of Problem |2j as inputs to 
Algorithm|2] From Lemma 5.7 we have E \\X' — Y\\ < J2i j 'oiij maXj^ecSajeCJ 11^^ ~ 2/11- Let A{Cr) be an optimal 



solution to Problem jlj over partition Cr, i.e. Ufir) — "^ij QLij ^^-y^c^x^c^ \\x — 2/11- Note because A{Cr) is the 
optimal solution to Problem |2j we also have E \\X' — Y\\ < 'J2ij n'ii'^^yec\xeci \\x — y\\ - Finally, we have that 



1(C) < W{ipp,,ipp) for any partition C (see proof of Lemma 5.2 in Appendix \B\. Combining these results, we 
obtain 



E\\X' -Y\\-W{^B,ifp) <J2q 

ij 



max \\x — y\\ — mm \\x — y\\ 
yec yec 
xec xec 
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Proof of Lemma 5.9 We first focus on the case d > 3. The proof rehes on the characterization of the length 
of the bipartite matching produced by Algorithm |3] (which also bounds the length of the optimal matching). By the 
triangle inequality, the length Ly[{Qn) of its matching is at most the sum of the matches between X and X' , plus 
the distances from the sites in y to their shadows in X' , i.e. 



LyiiQn) < Lyi{{X' , X)) + Lyx' 



(19) 



where Lyx' = TI^y x')eM W-^' ~ ^li- subtracting on both sides of equation ([T9| the term nW((pu, fp), and 
dividing by n^^^^'^, we obtain 



Lm{{X',X)) 



^nEWX' -Y\ 



n 



l-l/d 



O 



where the last equality follows from Lemma 5.8 Lemma 5.6 allows us to apply equation Q to Lyi{{X' , X)), and 
so the limit of the first term is 

Lu{{X',X)) 



lim 

n— >-+oo 



almost surely. We observe that riE||X' — y|| is the expectation of Lyx', and so the second term goes to zero 
almost surely (absolute differences law. Section II-Di. The resolution function of Algorithm [3] ensures that the third 
term vanishes. Collecting these results, we obtain the inequality in ([8]l with (f> = ipp. To complete the proof for the 
case d > 3, we observe that Algorithm |2] could be alternatively defined as follows: the points in X generate a set 
y' of shadow sites; the intermediate matching is now between y and y' . One can then prove results congruent 



with the results in Lemmas 5.6 5.7 and 5.8 By following the same line of reasoning, one can finally prove the 
inequality in ^ with (j> = (^d- This concludes the proof for the case d> 3. The proof for the case d = 2 follows 
the same logic and is omitted in the interest of brevity. ■ 
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